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Abstract 

A novel hybrid Monte Carlo transport scheme is demonstrated in a scene with solar 
illumination, scattering and absorbing 2D atmosphere, a textured reflecting mountain, 
and a small detector located in the sky (mounted on a satellite or a airplane). It 
uses a deterministic approximation of an adjoint transport solution to reduce variance, 
computed quickly by ignoring atmospheric interactions. This allows significant variance 
and computational cost reductions when the atmospheric scattering and absorption 
coefficient are small. When combined with an atmospheric photon-redirection scheme, 
significant variance reduction (equivalently acceleration) is achieved in the presence of 
atmospheric interactions. 
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1 Introduction 

1.1 Motivation and Background 

Forward and inverse linear transport models find applications in many areas of science 
including neutron transport [1, 2, 3], medical imaging and optical tomography [4, 5], 
radiative transfer in planetary atmospheres [6, 7, 8] and in oceans [9, 10], as well 
as the propagation of seismic waves in the solid Earth [11]. In this paper, we focus 
on the solution of the forward transport problem by the Monte Carlo (MC) method 
with, as our main application, remote sensing (an inverse transport problem) of the 
atmosphere/surface system [12]. In our demonstration, light is emitted from the Sun 
and propagates in a complex environment involving absorption and scattering in the 
atmosphere and reflection at the Earth's surface before (a tiny fraction of) it reaches 
a narrowband detector, typically mounted on a airplane or a satellite. 

The integro-differential transport equation (1) may be solved numerically in a va- 
riety of ways. Monte Carlo (MC) simulations model the propagation of individual 
photons along their path and are well adapted to the complicated geometries encoun- 
tered in remote sensing. Photons scatter and are absorbed with prescribed probability 
depending on the underlying mcdivim. The output from the simulation, e.g., the frac- 
tion of photons that hit a detector, is the expected value of a well-chosen random 
variable. These simulations are very easy to code, embarrassingly parallel to run, and 
suffer (in principle) no discretization error. The drawback is that they can be very slow 
to converge. MC methods converge at a rate (variance /NY/"^ where N is the number 
of simulations, and the variance is that of each photon fired. In remote sensing, the 
(relative) variance is high in large part because the detector is typically small and thus 
most photons are not recorded by the detector. In order to be effective, even in a 
forward simulation, MC methods must be accelerated. 
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One approach to speedup MC simulations is to use quasi-Monte Carlo methods, 
which steepen the convergence rate from ~ N~^/'^ to a more negative exponent. How- 
ever, most MC speedup efforts focus on reducing the variance of each photon. See 
[2, 3] or the review of more recent work on neutron transport in [13, 14, 15] and on 
3D atmospheric radiative transfer in [16, 17]. See also [18] for a thorough introduction 
to the MC techniques, including variance reduction, used in computer graphics. In 
problems with a small detector, this is achieved by directing photons toward that de- 
tector, and re-weighting to keep calculations unbiased. When survival-biasing is used, 
photons have their weight decreased rather than being absorbed [2, 3].^ Often, one 
uses some heuristic (such as proximity to the detector), or some function to measure 
the "importance" of each region of phase space. In splitting methods [2, 3], the photon 
is split into two or more photons upon identifying that a photon is in a region of high 
importance. The weight of each photon is then decreased proportionately. Propagating 
many photons with a low weight is not desirable, therefore splitting is often accompa- 
nied by Russian roulette. Here, if a photon enters a region of low enough importance, 
then the photon is terminated with a certain probability, i.e., high chance of absorp- 
tion if the weight is low; in the rarer alternative outcome of the Bernoulli trial, the 
weight is increased to keep the simulation numerically unbiased. So there is typically a 
slight cost in variance to improve efficiency (by terminating low- weighted trajectories). 
Typically a weight window is used to enforce regions of low/high importance. Source 
biasing techniques change the source distribution in order to more effectively reach the 
detector. More generally, the absorption and scattering properties at any point can be 
modified, provided photons are re-weighted correctly. 

It has long been recognized that the adjoint transport solution is a natural impor- 
tance function [19, 2, 3, 20, 21, 22, 14, 23, 15]. One can use approximations of the 
adjoint solution — typically a coarse deterministic solution — to reduce variance. The 
result is a hybrid method (deterministic k, MC). The AVATAR method uses an adjoint 
approximation to determine weight windows [22]. The CADIS scheme in [14] uses an 
adjoint approximation in both source biasing and weight-window determination. An 
adaptive technique that successively refines the solution in "important" regions, using 
the adjoint to designate such regions, is described in [24, 25]. In [19, 2, 3, 20], a zero- 
variance technique is outlined that uses the true adjoint solution to launch photons 
that all reach the detector with the same weight ... which happens to be the correct 
answer. This method is of course impractical since determining the exact adjoint so- 
lution everywhere is harder than determining some specific integral of that solution, 
which is usually the goal of a MC simulation. The LIFT method [20, 21] therefore uses 
an approximation of the adjoint solution to approximate this zero-variance method. 

We adapt the zero-variance technique to the particular problem we have at hand; 
see Fig. 1 for the type of geometry considered in this paper. The problem we consider 
has a fixed, partially-reflective, complex-shaped lower boundary, and relatively large 
mean-free-path (MFP) in the sense that a large fraction of the photons reaching the 
detector have not scattered inside the (optically thin) atmosphere. Calculation of the 
approximate adjoint solution used to emulate zero- variance techniques is difficult and 

-'^Notc the somewhat confusing terminology: On the one hand, a method is statistically biased if the 
expected outcome is not the intended one. On the other, the practice of re-directing photons in favorable 
directions and/or reducing the number of scattering events is also called biasing. In the latter case the 
photon has its weight adjusted so that the simulation is unbiased. 
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potentially very costly. What we demonstrate in this paper is that partial, "localized" 
(in an appropriate sense) knowledge of the adjoint solution still offers very significant 
variance reductions. More specifically, we calculate adjoint solutions that accurately 
account for the presence of the boundary but do not account for atmospheric scattering 
(infinite MFP limit). The computation of the adjoint solution thus becomes a radiosity 
problem with much reduced dimensionality compared to the full transport problem. 
This, of course, can only reduce variance in proportion to the number of "ballistic" 
photons that never interact with the atmosphere. When combined with simple rules 
for allowing atmospheric scattering and sending some photons directly from the atmo- 
sphere to the detector, our hybrid method yields very significant variance reduction 
at relatively minimal cost. Furthermore, the methodology studied is applicable when- 
ever any method is available to deterministically pre-calculate flux over any subset of 
paths. For instance, complex propagation of light in clouds and its importance could 
be pre-calculated locally and incorporated into the MC simulations in a similar fash- 
ion. This "modular" approach to the description of the adjoint solution is well-adapted 
to the geometries of interest in remote sensing and avoids complicated, global (hence 
expensive) deterministic calculations of adjoint transport solutions. Our treatment of 
the reflecting boundary described in detail in this paper is a first step toward modu- 
lar adjoint transport calculations and their variance reduction capabilities in remote 
sensing. 



Figure 1: Mountain (1 — cos^ x shape), cloud, sky, and detector. Dot size indicates relative 
adjoint flux strength. Large dots on right-hand-side are the detector (dot size is down-scaled 
for detector). Dot size on mountain indicates that portions of the mountain are shaded from 
the detector, and that the surface albedo is varying. See section 3.1 for specifics, as used in 
the present study. 

The rest of the paper is organized as follows. In sections 1.2 and 1.3 respectively, the 
physical problem and statistical formulation are described. In section 1.4, the analog 
and survival-biased MC algorithms are presented. In section 2, the surface adjoint 
importance (SAI) and regularized SAI methods are introduced. These are the hybrid 
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adjoint-based methods at the core of this work. In section 3, numerical estimations of 

variance reductions and computational spcedups are given. For an expanded exposition 
of techniques and analyses from a mathematical standpoint, we refer the interested 
reader to [26]. Finally, we summarize our findings in section 4 and conclude with 
thoughts about potential applications in remote sensing science. 



1.2 Problem Setup 

Our setup is photon transport in a domain R cW^ {d = 2, 3, with d = 2 in our present 
demonstration) described in Fig. 1. The outward normal to the domain boundary dR 
at position r is denoted Vr- R is the atmosphere, the sky/mountain/sides/detector 
constitute dR. We have one small detector located on the right-side of dR, and the 
goal of simulations is to estimate the photon flux through the detector. Since photons 
are monokinetic, propagation direction w is a unit vector in the sphere S*^"^ embedded 
in d-dimensional space (unit circle for d = 2). 

We model radiance (a.k.a. specific intensity or angular photon flux density) /(r, v) 
in our medium with a boundary source distribution Q. I obeys the following integro- 
differential transport equation and boundary condition: 

V ■ V/(r, v) + a{r)I{r, v) = KI{r, v) 

/(r,.) = 4^ + ^, redR,.ndvUr<0, 

\l'r ■ v\ \Ur ■ v\ 

the integral operators being defined by kernels 

Kf{f->v) = <7s(r) / p{r,v'^v)f{r,v')dv', r&R 

Kf{r,v) = a(r) / P{r,v'^v)\i'r ■ v'\f{r,v') dv' r E dR, and t; ■ t'r < 0. 

Jl/r-v'>0 

(2) 

The extinction coefficient, a.k.a. total cross section (per unit of volume) C7(r) is 
the sum of the intrinsic absorption coefficient/cross-section aa{r) and the scattering 
coefficient/cross-section o"s(r). For the partially reflecting boundary condition (viewed 
here as a surface scattering), a{r) is the local value of the albedo. Both volume 
{p{r,v'^v)) and surface {P{r,v'—^v)) phase functions are normalized (J pdv' = 1). 
Since the transport problem is linear, we use a normalized boundary source, i.e., 

/ / Q{r,v)dn{r)dv = l (3) 

JdR Jur V<0 

where d/i(r) is the appropriate measure on the {d — l)-dimensional boundary. 

Our detector measures photon flux and is described by a "response function," 
g{r,v)\i'r ■ v\, where g{r,v) is zero everywhere except when r is in the physical de- 
tector (its aperture or "pupil") and v points out of the boundary. Where ^ it is 
constant and, furthermore, it is normalized so that J g{r, v) dr dv = 1. The goal of our 
Monte Carlo method is to compute the detector's signal 



/ / 

J dR J Vr-V 



g{r, v)\i'r ■ v\I{r, v) d/u(r) dv. (4) 

>o 
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In the present study, any v can contribute to the radiometric signal measured at r. To 
model an imaging detector, direction space would be limited to a finite field-of-view 
that would in turn be subdivided into individual "pixels." 
For future use, we define the function 

r f\r-r'\ ] 

Ea{r,r') ■. = eyip<- J a{r + tr' - r)dt> , 

where r' — r := (r' — r)/\r' — r\. Physically, it describes the probability of direct 
transmission of light from point r to point r' (or vice- versa), that is, without suffering 
any collision. 

1.3 Statistical formulation/notation 

The measurement defined formally in (4) is approximated in a Monte Carlo simulation 
by estimating an average 

1 ^ 

n=l 

where the a;„ are photon paths (w = (ro, ri, . . . , rfc)) generated by the "analog" chain 
(meaning analogously to real-life photon travel, cf. Algorithm 1), and the relevant 
indicator function is l£)(w) = 1 if the path hits the detector (hence the subscript D), 
and = otherwise. The paths are random variables and, with E {•} denoting statistical 
expectation, we have 

E{Sn} = E{1d} = 'P[D], 

where the notation P[D] emphasizes that this is a probability of hitting the detector. 
We also have P[D] equal to the desired measurement or signal in (4). For finite N, Sn 
is not equal to P[D] exactly. The mismatch is quantified in a statistical sense through 
the variance 

Var {Sr,} : = E {(5^r - P[D]f} = ^ = 

since all of the events ojn contributing to the Sn estimator are independently drawn. 

Rather than Sn, one may generate paths according to some modification of real-life 
photon travel and then estimate 



1 ^ dP^ 
P[D] ~ Tjv := - 5^ lB(a;„)-=-(a;„), 



n=l 

where the ratio dP^/ dP(c<;) is the ratio of the probability density of u in the analog 
chain to that in the modified chain. This importance sampling technique is widely used 
in statistics since often times Tn will have lower variance than Sn- Indeed, most of the 
variance reduction techniques mentioned in the introduction are of this type. In our 
algorithms we compute this ratio step- by-step and refer to it as a weight (modifier). 
So, rather than counting photons, we count weighted photons. 



6 



For future use we define the following (standard) statistical notations and conven- 
tion. First, we write u ^ Z//[0, 1] to indicate that n is a random variable uniformly 
distributed on the interval [0,1]. Second, a probability density such as 7r(x) can be 
denoted explicitly (e.g., 7r(x) = (27r)~^/^ exp {— in the case of the normal distri- 
bution with zero mean and unit variance), or it can be given up to a constant (since it 
must integrate to one). In this last case, we would write it{x) oc exp {— a;^/2}. 



1.4 Standard Algorithms 

We present here two basic algorithms for Monte Carlo transport. These are well known 
but we do this in order to demonstrate our notation. Algorithm 1 is often referred to 
as analog since the photons follow a path analogous to photons in the real world. 



Algorithm 1 Analog Monte Carlo Transport 
1: Choose a starting position/direction (ro, Vq) according to the sun's source density Q{r, v) 
2: Draw u ~ W[0, 1] and cast the photon along the ray ro + tvQ until E^^tq + tvo) < u. Call 

this point ri. If this does not happen before dR is reached then set ri to the boundary 

point at the intersection with the ray. 
3: if ri e -R then 

4: With probability (Js(ri)/cr(ri), the photon is not absorbed, and we select Vi using the 
probability density 



otherwise the chain is stopped. 
5: else if ri e dR then 

6: With probability a{ri) the photon is not absorbed, and we select vi using the proba- 
bility density 



otherwise the photon is absorbed and we stop the chain. 
7: end if 

8: Continue alternating casts and direction changes until either the photon is absorbed, 
escapes through the upper boundary ( "sky+sides" ) , or the detector is reached. 



For use in Algorithm 5 further on, we will need to know the probability density of 
the analog chain producing a path oj. This is given by 



vi ^ P{ri,Vo-?-Vi); 



D analog {ro, n) = QirQ, Vq) E^{ro, n) , 
Danaiogiro,ri,r2) = I^ana^ (^0 , ri)K"'"''°^ (h , r7^i)E^{ri,r2) 



and so on. Above is given by 




n e dR. 
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Algorithm 2 uses a trick known as survival-biasing since photons will survive (al- 
most) any interaction with the media. We do this by casting photons while ignoring 
intrinsic absorption. So, e.g., if a patch of media has aa large and < (Js/cr <C 1 the 
photon will almost never scatter there. Our weight is then Efj/E^r^. When the photon 
interacts with the surface, then so long as a > 0, we do not absorb but multiply the pho- 
ton weight by a. Another, slightly different but also common, survival-biasing method 
would cast photons in the same manner as analog, but would eliminate absorption and 
re-weight by as/fJ. So, e.g., if a patch of media has Ga large and < cts/ct ^ 1 the 
photon would likely interact with the media and scatter but not be absorbed there; its 
weight however would be reduced by a factor of (Ts/ct (known in the radiative transfer 
literature as the "albedo for single scattering"). 



Algorithm 2 Survival-Biased Monte Carlo Transport 
1: Choose a starting position/direction (ro,vo) according to the source density Q{r,v) 
2: Draw u ~ W[0, 1] and cast the photon along the ray rg + o until [tq + Ivq) < u. Call 
this point ri. If this does not happen before OR is reached then ri is the boundary point 
we have reached. Since we paid no attention to intrinsic absorption during the cast, the 
photon picks up a weight equal to 

K„(ro,ri) = — 

(ro,ri) 

3: if ri G -R then 

4: Select vi using the probability density 

5: else if ri G dR and a{ri) > then 

6: Select vi using the probability density 

vi ^ P{ri,VQ-^vi). 

Since we had no chance of boundary absorption, the photon's weight is multiplied by 
a{ri). 

7: else if ri G dR and a{ri) = then 

8: The photon is absorbed and we stop the chain. 

9: end if 

10: Continue alternating casts and direction changes until either the photon is absorbed, 
escapes, or reaches the detector. 



2 The Surface Adjoint Importance (SAI) Method 

The SAI method uses an approximation to the surface reflection problem to reduce 

variance coming from surface interactions. It ignores atmospheric effects and therefore, 
by itself, is statistically biased. In section 2.2 we pair it with other methods to produce 
an unbiased estimate of the detected flux. 
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2.1 Pure SAI 



Here we ignore atmospheric effects and demonstrate and develop a Monte Carlo metfiod 
tfiat sends photons from surface point to surface point and then to the detector. If at- 
mospheric effects are not present, and our deterministic solution was perfectly accurate, 
this method would have zero variance. 

The adjoint solution to transport may be developed by considering the L? adjoint 
of the integral solution to transport and reversing the role of the source and detector. 
Let P be the adjoint solution when only surface effects are present. We therefore have 

P{r, v) = a{r) / P(r, v^v')P{r+{r, v'),v') dv' + g{r, v). (5) 

Jl/rV'<0 

This adjoint solution corresponds (in a Monte Carlo viewpoint) to sending photons 
that start at the detector and travel backwards. Therefore, it will have its maximum 
at the detector. It will be higher in places that have a clear path to the detector. 
will be zero at places from which a photon cannot reach the detector. Our numerical 
solution of (5) is described in the Appendix. 

The pure SAI chain is defined by the steps described in the following Algorithm 3. 

Algorithm 3 Pure SAI 
1: Choose a starting position/direction (ro,fo) according to the modified source density 

Q'''%r,v) oc Q{r,v)r{r+{r,v),v). 

The photon picks up a weight Q{ro,Vo)/Q^"-^{ro,Vo) 
2: Cast the photon until it hits the opposing boundary at point/direction {ri,Vi) — 
(r+(ro, I'l), I'l). The weight is multiplied by 

1 

3: Change direction according to the density 

where Csai(r) is a normalization factor depending only on r e dR. Since we did not 
account for boundary absorption, the photon weight is multiplied by a{ri). The modified 
direction change must also be taken into account and therefore, in addition, the weight 
is multiplied by 

4: Cast the photon until it hits the opposing boundary. If it hits the detector, stop and 
record a hit. Else, repeat step 3. 



For use in Algorithm 5, we will need to compute the probability density of a path 
Lo = (ro, . . . ,r^) being generated by pure SAI. This is simply the denominator in the 
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corresponding weights. Denoting this by Dg^i we have 

Dsai{^) = 0, if rj € R for any j, 
and for paths such that rj G dR for all j, we define Dgai recursively (with Vj := 

Dgai{ro,n) = S^^'iro^vo) 
Dsai{ro,ri,r2) = Dsai{ro,ri)K^'''-{ri,vo -^vi), 
Dsai{ro, ■■■,rk) = Dsai{ro, ■ ■ . ,rk-i)K^''\rk-i,Vk-2 Vk-i). 

Remark 2.1. 

• A discretized version of the density Q*"' is pre-computed using the (discrete) 
solution see section 3. This means that we can pre-compute the normalization 
factor Csai- The discrete density Q^^^ will be defined at a number of points (r^, vq) 
where vq is the anti-solar direction. We use the density to decide on a center point 
Tj, and then perturb the starting point by a small (random) amount to eliminate 
discretization effects in the final solution. 

• The direction change pdf is also pre-computed and stored as a discrete pdf over 
angles. We use the pdf to pick a direction center vj and then perturb to obtain 
the new direction. 

• Up to numerical error one can see that Csai{r) = a{r). Indeed, dividing (5) 
through by P (r, v) we have 

I = air) I P(r,v^v') ,V C + \ - 

So away from the detector g{r,v) = and the integral is therefore equal to 1. 

• That this method is biased is easy to see: If a region of the atmosphere has non- 
zero scattering, then it would be possible (in the analog world) to scatter from 
that point to the detector. This type of interaction is not allowed in a pure SAl 
world. 

This is an implementation of the zero variance adjoint-based chains studied in 
[2, 3, 20] in the special case where atmospheric effects are not present. Hence (dis- 
regarding numerical error), this would be a zero- variance method were atmospheric 
absorption/scattering absent. 

We verify this claim numerically by testing the method in simulations without 
atmospheric effects. See Fig. 2 where this is tested with both a flat terrain and a 
curved "cos^" mountain. The curved mountain increases variance since discretization 
does not allow the function r+{x,v) to be implemented perfectly. 

2.2 Regularized SAI 

Here wc use the SAI chain as part of a larger unbiased chain. Since Algorithm 3 does 
not generate paths following all possible interactions, wc must supplement it with an 
algorithm that does. We then use a number qg G [0, 1] to determine the fraction of 
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Figure 2: Left: When a flat mountain is used, variance ~ 0{h^-^) where h is the discretization 
parameter. Right: On the curved boundary discretization effects are more prevalent and 
convergence is slower. 



photons that travel according to Algorithm 3 (this fraction = I — Qs), and what fraction 
according to the supplemental algorithm. 

Before describing the regularized SAT algorithm, we present the supplemental Al- 
gorithm 4 dubbed "heuristic scattering adjustment." It is a survival-biased algorithm 
in the sense that no absorption occurs within the atmosphere, or at boundary points 
(unless the boundary point had a = 0, e.g. the sides/sky). It also makes use of a 
simple scheme to direct a fraction of atmospheric interactions toward the detector. No 
claim is made to the optimality of this re-direction (it is similar to the technique of 
local estimation [27, 16]). We use Algorithm 4 since it is simple to understand and 
illustrates the dramatic decrease in variance that can be achieved when two methods 
(SAI and heuristic) are used together in Algorithm 5 (see also section 3). 
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Algorithm 4 Heuristic Scattering Adjustment with Parameter e [0, 1 



1: Choose a starting position/direction according to the standard source density Q{r,v) 
2: Cast the photon as in Algorithm 2 until it hits the opposing boundary or interacts with 

the atmosphere at ri. The photon picks up a weight equal to E„{rQ,ri)/ Efj^{rQ,ri) 
3: if ri e i? then 

4: With r^o = "the midpoint of the detector" , compute 

qheu[ri,Vo) :^l-{l-qy)- 



\\p{ri,vi^-)\\L'^ 

With probability 1 — g/ie« draw Vi from a uniform distribution of directions pointed 
toward the detector (we call this fviri,Vi)), and with probability Qheu draw Vi from 
p{ri,Vo-^-). The weight is multiplied by 

as{ri)p{ri,Vo-^Vi) 

if ri e K. 
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(1 - qheu)fv{ri,Vi) + qheuP{ri, Vq-^Vi) ' 

5: else if ri e dR and a{ri) > then 

6: pick a new direction according to the density P{ri,Vo-^Vi). The weight is multiplied 

by a(ri). 

7: else if If ri G 9/2 and a(ri) = then 
the photon is absorbed and we stop, 
end if 

Continue in tliis niainier until al)S()r])ti()ii or tlie dotoctor is reacried 



So at every scattering event, the weight is modified by a ratio of either a{r)P or 
(Js{r)p to K^^'^ where 



(1 - qheu)fv{r, v) + qheuP{r, v^'), r e R 

P{r,v^v'), redR. 



For use in Algorithm 5 we will need to compute the probability density of a given 
path generated by Algorithm 4. This is simply the denominator in the corresponding 
weight. Denote this by -D/iew('"0! ''ij • • • i '^fc); which we define recursively by (with Vj := 

Dheu{ro,ri) = Qiro,vo)Ea,{ro,ri), 
Dheu{ro,ri,r2) = Dheu{ra,ri)K^^''{ri,VQ vi)E^^{ri,r2), (6) 
Dheuiro, ...,rk) = Dheuiro, ■ ■ . ,rk-i)K'^'''^{rk-i,Vk-2 -> Vk-i)Eas{rk-i,rk), 

and so on. 

We now present Algorithm 5, the regularized SAI algorithm that combines pure 
SAI (Algorithm 3) with the heuristic scattering adjustment (Algorithm 4). Note that 
any unbiased algorithm may be combined with pure SAI in a similar manner. 
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Algorithm 5 Regularized SAI with parameters qs,qv £ [0, 1] 
1: With probability 1 — qs, generate a path according to Algorithm 3. With probability q, 

generate it according to Algorithm 4. 
2: The weight of the path cu = (ri, . . . , r,-) is 

D analog ('^) 

(1 - qs)Dsai{uj) + qsDheu{(^) ' 



Algorithm 5 uses SAI to produce paths that interact only with the surface. One 

could easily devise other algorithms that send paths via the heuristic chain, and once 
paths interact with the surface they use the SAI chain. This could reduce variance 
further, but we choose not to study this in order to simplify the presentation. 



3 Numerical Results 

3.1 Parameter choices in numerical simulations 

In the assumed d = 2 transport space, we have r = (x, y), where x increases from 
left to right in Fig. 1 and y increases from bottom to top; r = (0, 0) is the point at 
the bottom of the valley. For directions, we have v = v((f>) = (cos ^, sin 0) where (f) 
increases counterclockwise from the x > axis. 

In the simulations performed with a = (no atmospheric interactions), we used 
both a flat surface (so that our domain was [— tt, tt] X [2, 4]) and a "cos^" surface (Fig. 1). 
We swept h, with 0.002 < h < 0.2. Wc did not use any heuristic scattering adjustment 
(g„ = 1.0). In all cases, we assume an isotropic (Lambertian) redistribution by diffuse 
surface reflection. This leads to the following surface scattering phase function and 
assumed surface albedo distribution: 

P(r v^v') (X i.^^^ ' t'r-'y'<0_ ^c^\ _ ff; |x| < 2.5 
\ 0, otherwise ' \0, otherwise ' 

The cutoff I a; I < 2.5 was done to simplify the coding (allowed us to use one simple 
routine for all values of h), and has no theoretical consequence. The source was mono- 
directional (p = —T^/^, and given by 

In the simulations involving atmospheric interactions (a > 0), we used a cos^ type 
surface. We compute speedup in a variety of cases. The mean-free-path MFP= 
was varied as well as Qs, h, and g^. We swept 0.002 < h < 0.15. In all cases the 
atmospheric scattering coefficients were constant with as = 2aa (hence (Js/cr = 2/3). 
The atmospheric scattering was given by 

p{r, v^v') ca 1 + {v ■ v')'^, 

which mimics a molecular (Rayleigh) in d = 2. The other coefficients were chosen to 
have features (in this case oscillations) on a scale coarser than the fine values of h, and 
finer than the coarse values. 
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The surface albedo was chosen to be quite complex (significantly different than the 
flat surface/constant reflection commonly used). The phase function P was as before 
(Lambertian) , but a is given by 

( \x\> 2.5, 

a{r) = I 0.75 + 0.25 sin(27rx/0.05) 1< x < 2.5, 
1^0.35 + 0.25sin(27rx/0.05) -2.5 < x < 1, 

using the same inconsequential cutoff \x\ < 2.5. Off the mountain there was no scat- 
tering (perfectly absorbing boundary). 

The source was mono-directional = — 7r/2 and given by 



Q{r,v{-TT/2)) oc 



1 + 0.25 sin(27rx/0.07) \x\ < 2.5, y = 4 
otherwise 



3.2 Speedup (figure of merit) 

We start by defining our figure of merit used to compare the different algorithms. 
We take the viewpoint that each algorithm produces a sequence of paths and 
corresponding random variables ^(w"^) equal to the product of Id^lo^) times the weight 
that the photon picked up along the way. To distinguish different methods we write 
for analog, ^gb for survival-biasing, ^gai for pure SAI, ^^eu for heuristic scattering, and 
for the regularized SAI method. 

For all of these methods, define the approximation after N random draws 

1 ^ 

n=l 

For ^ equal to any of the above methods, In{0 unbiased estimator of E{^} = 

P[D], i.e., the probability of a detector hit. 
The RMS estimation error e is given by 

e{o ■■ = v^{\iN{o-nDW} = ]m^- 

For a given error level e, the required number of MC draws is then N(e, ^) := Var {^} /e. 
The required simulation time T{£,^) for one estimation of P[D] is given by 

T{s, : = ToiO + rm = MO + 1^^)^^, 

where ro(^) is the time needed to compute the deterministic adjoint solution (e.g. at 
level h when ^ = ^/j), and r(^) is the expected time for one draw using the appropriate 

measure for the random variable ^. We foresee the use of SAI in situations where the 
boundary remains fixed, but the atmosphere changes (due to, e.g., moving clouds over 
a fixed surface). We therefore consider the time for m simulations using one boundary, 

T{s, m) : = ToiO + mriON = MO + m^^^^^^. 



14 



Schemes may be compared with the ratio 



T{s,^i,m) _ £^ro(gi) + mr(6)Var{a} 
Tie, 6, m) s^To{^2) + mT(6)Var {6} ' 

For a deterministic approximation of we expect ro(^) C{$,)h~^^'''~^\ We in 
fact measure (with d = 2) To{^h) ~ 0.017h~^. Our "benchmark" scheme is survival- 
biasing. Since ^sb requires no deterministic solution, the relevant ratio (and our figure 
of merit) is 

Speedup(^5,£:,m) • — 



(f)^C + mr(e,)Var{CJ 



We measured speedup when either m = 10 or, formally, m = oo ("Ignoring deter- 
ministic solve"). 



3.3 Variance reduction 

Here we analyze the variance of the SAI chain in the presence of atmospheric interac- 
tions. Note that even when the error \P[D] — {P, S) \ is high, we still get good variance 
reduction. See Fig. 3. This emphasizes the point that the quality of the deterministic 
solve is not so important in a modular scheme. 




Figure 3: \P[D] — (J*, S')|/P[D] is generally lower for smaller h. However, speedup is still 
very good even for large h. Diam is the maximal diameter of the simulation domain R 



Our implementation swept both qs and qy. As expected, we see decreasing speedup 
with increasing atmospheric scattering strength a. See Fig. 4. 

It is important to note that use of adjoint-enhanced surface scattering, and heuristic 
atmospheric scattering (g^ < 1, g^, < 1) together is especially helpful. In fact, even with 
a small MFP = 1.3-Diam (Diam is the maximal diameter of the simulation domain R), 
we realize good speedup when qg = 0.9, q^ < 1. Note that if either qg = 1 or qy — 1 (so 
no use of either SAI or heuristic scattering adjustment), speedup almost disappears. 
This is slightly counter-intuitive but may be explained as follows: Each method (SAI or 
heuristic) significantly increases the number of paths in two significant classes (surface- 
only and atmospherc-to-dctcctor) . Therefore, variance from these path-classes is all 
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cos^ mountain, MFP = 16xDiam cos^ mountain, IMFP = SxDiam 




Figure 4: Speedup when using both surface adjoint approximation P (with parameter 
and heuristic atmospheric scattering (with parameter q^) 



but eliminated. Supposing each of these path-classes accounts for 2/5 of the total 

paths reaching the detector, by themselves they can only reduce variance by a factor 
of 1/(1 — 2/5) = 5/3. However, together they can reduce variance by a factor of 
1/(1 -4/5) = 5. 

As one can see, selection of the parameters Qs and makes a significant difference 
in the resultant variance. We provide some heuristics here and refer the reader to [26] 
for more details. When — >■ most of the photons will travel on the surface only. 
The photons that take a route prescribed by the heuristic chain must then carry an 
additional weight = I/qs to compensate for this. For this reason, picking Qs too small 
results in increased variance. A similar argument holds for q^. That the optimal qg is 
so close to 1 (and greater than the optimal q^) can also be explained by the fact that 
paths interacting exclusively with the surface are less likely to occur (in the analog 
world) than those interacting with the surface and atmosphere. 

4 Conclusion and Outlook 

A novel method for Monte Carlo transport was presented that uses an approxima- 
tion of the adjoint (ignoring atmospheric effects) to reduce variance in simulations, 
equivalently, accelerate convergence to a specified accuracy. This algorithm, the Sur- 
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face Adjoint Importance (SAI) method, may be combined with any unbiased method 
to significantly reduce variance coming from surface interactions when the overlaying 
atmosphere is optically thin. If it is combined with a method that reduces variance 
coming from atmospheric interactions, significant overall variance reduction is achieved. 
The implementation is relatively simple, requiring only an approximate adjoint trans- 
port solver for the boundary which adds virtually no overhead to the Monte Carlo 
computation time. 

A possible application of this kind of accelerated Monte Carlo modeling in remote 
sensing is to address "adjacency" effects caused by highly variable terrain, including 
built environments (urban canyons). The standard adjacency effect is observed when an 
aerosol layer of moderate optical thickness mixes in an imaging detector's pixel light 
that has been reflected off surface elements with contrasting albedos in neighboring 
pixels. This is now a solved problem in the case of a variable-but-flat surface under a 
uniform atmosphere [28]. However, adjacency effects caused by non-flat terrain arc only 
beginning to be explored, particularly in the thermal IR (where Q{r, v) is determined 
by temperatures and emissivities). 

On a broader scale, our work is an illustration of a modular approach to variance 
reduction whereby different interactions are handled separately and then pieced to- 
gether in an unbiased manner. Specifically, these different interactions could be pieced 
together as in Algorithm 5. 

For instance one can envision a "cloud" module where radiation transport inside the 
cloud (dominated by multiple scattering) is treated off-line in some judicious approx- 
imation, and then incorporated into complex scene simulation. In applications driven 
by surface property retrievals from remote sensing data, efficient modularized Monte 
Carlo modeling would open the door to advanced atmospheric compensation schemes 
with broken-cloud capability. This is another wide open frontier recently explored in 
[29]. 
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A Appendix: Numerical solution to the adjoint 
problem 

Here, at discretization level h, we approximate P . 

To simplify computation of our numerical solution we make the assumption 

P{r,v^v') = l^^.y>Q{r,v)K{r,v'), 

and recall that g{r,v) = goir) = constant whenever Ur ■ v > so that g{r,v) = goir)- 
The result is that is then a function of position only. This significantly improves the 
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speed of solving the adjoint problem, as well as the memory requirements for using it. 
Theoretical results in this paper do not need this assumption, which we make here as 
a matter of convenience. 

We will now discretize the coefHcients and approximate the integral operator ap- 
pearing on the right hand side of (5), denoted now by T. For r\ G dR, 

Tr{ri,vi) = a{ri) / K{ri,V2)P{r+{ri,V2),V2)dv2. 



■V2<0 

Notice that Tf is function depending only on r, and in fact only on the boundary 
values of /. Since g depends only on r, = J2'h=o'^''S ^^^^ depend only on r and 
whether or not Uj. -v > 0. We thus define 

(p{r) : = P{r,v), r G dR, Ur ■ v > 0. 

We find that ip : dR — )• R satisfies the equation 

(p = A(p + go, Af{ri):=a{ri) K{ri,V2)f{r+{ri,V2))dv2. 

J Uri ■V2<0 

In discretizing this operator, and integrals over directions in general, we use the 
change of variables, 



/ 



f{r+{r, v), v)dv= / /(r', v)d^N{r, r') d/i(r'), 

I JdR 



yr-v<0 JdR /• -■ \ 

The term d^N is normal derivative (at r) of the free-space Green's function for the 
Laplacian. One can show (sec, e.g., the section on double-layer potentials in [30]) that 
for r,r' G dR, Vr ' {f' — r) ^ \r' — r|^. Therefore it is in fact an integrable function. 
When d = 2 it is moreover bounded. 

We now discretize the operator A. First split the boundary into non-overlapping 
segments {dRj}^^^ with dRj centered at Vj, with length \dRj\ < h. Denote by Rf the 
(orthogonal) projection of / onto the space of piecewise constant functions (constant 
on each segment dRj). We also think of Rf as a vector in and Rfj its components. 
Then, after the change of variables (A.l) we have (at gridpoint r,) 

Af{ri) = a{ri) / K{ri,f'^)dt,N{ri,r)f{r)diJ,{r) 
JdR 

~«(n) \9Rj\K{ri,f^^i)d^N{ri,rj)f{rj) 

0<j<Np-l i-^-'^) 

j 

This implicitly defines the matrix A^. 

We now define our discrete approximation to (p as the piecewise constant function 
(vector) tp^ solving 

(^'^ = A'^ip'' + Rg. (A.3) 
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We then define approximations !^ P, 

/^(r, v) : = (p''{r), r e dR, Ur-v>0. (A.4) 

Note that, in our implementation, we have chosen to represent angular integrals as 
integrals over the boundary. This works for two reasons. First, as our adjoint solution 
depends only on position it is convenient to evaluate these sums. Second, if instead 
a discretization were chosen that was uniform in angle, then (with only finitely many 
angles) one would often miss the (small) detector in evaluation of the integral. 
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